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Preface 

On  Attgnit  28,  1978,  the  fint  artificial  satdlite  was  launched  into  an  unstable 
periodic  orbit  about  the  interior  Ubration  point  (or  Lagrange  point)  between  the 
sun  and  the  earth.  It  was  called  International  Sun-Earth  Explorer,  ISEE-3,  and 
its  purpose  was  to  make  accurate  measurements  of  the  magnetic  phenomena  in  the 
solar  wind.  Proposed  by  Fbrqnhar  (2),  it  was  onr  first  Ubration  point  satdUte, 
demonstrating  what  Floqnet,  Poincari  and  Lagrange  could  only  think,  write  and 
dream  about. 

Other  missions  using  similar  orbits  have  been  proposed,  and  the  purpose  of 
this  thesis  is  to  discuss  the  control  of  satelUtes  in  such  an  orbit  Floquet  theory. 
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Abstract 

A  periodic  “halo’'  orbit  which  eziiti  about  the  interior  Lagrange  point  for  the 
Earth-ran  system  was  decomposed  using  Floqnet  theory  into  modal  irariaUes,  which 
are  dynamically  decoupled  subspaces  for  the  six  degree  of  freedom  system.  Modal 
control  consisted  of  evaluating  the  diverging  mode  and  maneuvering  to  counteract 
its  divergence.  In  the  unperturbed  system  this  was  successful.  Control  costs  were 
low,  and  the  significance  is  that  the  controller  did  nothing  to  suppress  modes  that 
were  oscillatory  or  converging.  First  order  perturbations  by  the  moon  and  Earth’s 
eccentricity  were  examined.  The  effect  of  the  moon’s  motion  allowed  the  scheme  to 
operate  with  reasonable  control  costs,  but  the  effect  of  eccentricity  caused  divergence 
in  spite  of  the  controller. 


Modal  Control  of  a  Satellite  in  an  Unstable  Periodic  Orbit  About  the 
Earth-Sun  Interior  Lagrange  Point 


1.  Introduction 

It  ii  desirable  to  have  a  satellite  upstream  of  the  Earth  to  study  the  solar  wind 
and  to  warn  us  of  variations  in  solar  activity.  The  mathematics  of  the  restricted 
three-body  problem  predict  an  equilibrium  point  in  the  space  between  Earth  and  the 
sun.  Ideally,  a  satellite  placed  at  this  point  could  literally  remain  there  indefinitely. 
However,  a  satellite  at  the  interior  Lagrange  pdnt  (or  Xa)  is  directly  between  earth 
and  sun,  so  from  any  observation  station  on  Earth,  it  would  appear  in  the  center 
of  the  sun’s  disk.  This  is  a  poor  place  from  which  to  receive  the  very  low-power 
transmissions  necessary  for  data  reception  and  tracking.  We  can,  however,  orbit  a 
satellite  about  the  Lagrange  point,  as  was  successfully  demonstrated  in  1978  with 
ISEE-3  (6).  If  a  simple  periodic  orbit  can  be  found,  it  is  called  a  Halo  orbit.  is 
an  unstable  libration  point  in  the  three-body  problem  (8:  271)  and  the  halo  orbit  in 
the  vicinity  of  X3  that  we  used  for  this  study  is  also  unstable.  Any  satellite  placed 
into  such  an  orbit  would  require  station-keeping.  Without  some  control  input  over 
time,  a  satellite  placed  near  such  an  orbit  will  not  remain  near,  but  quickly  wander 
off  into  space  surprisingly  far  from  the  intended  path.  Even  if  we  could  iigect  a 
satellite  exactly  into  a  halo  orbit  and  not  just  near  one,  it  would  not  remain  there 
indefinitely  since  the  universe  we  live  in  is  not  completely  described  by  the  restricted 
three-body  problem.  This  is  analogous  to  trying  to  balance  a  pencil  indefinitely 
on  its  point.  Even  if  you  placed  it  perfectly  in  balance  and  there  were  no  wind  or 
vibrations,  eventually  some  ridiculously  small  force  such  as  the  coriolis  effect  or  the 
moon’s  gravity  would  knock  it  over.  This  paper  will  concentrate  on  the  instabilities 
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inlieient  in  tke  orbit  and  tke  perturbing  effiecti  of  the  Earth’s  eccentricity  and  the 
moon. 

The  station-keeping  of  ISEE-3  was  accomplished  by  calculation  of  a  nominal 
periodic  orbit,  a  reference  orbit,  from  which  it  was  allowed  to  stray  only  a  certain 
amount  before  its  small  thrusters  were  fired  to  kick  it  back  closer  to  the  reference 
orbit.  There  may  be  a  better  way.  In  this  paper  we  will  use  the  theory  of  Floquet 
to  demonstrate  that  while  the  orbit  is  unstable,  only  one  of  its  six  Floquet  modes  is 
imstable.  It  is  our  goal  to  suppress  only  this  mode,  while  leaving  the  stable  modes 
and  the  oscillatory  modes  to  their  own  devices.  The  stable  modes  will  tend  to  bring 
the  satellite  closer  to  the  reference  orbit  in  the  long  run,  even  if  the  orbit  itsdf 
appears  to  be  diverging  at  any  given  time,  and  likewise,  the  oscillatory  modes  will 
not  diverge,  but  wobble  back  and  forth  across  the  nominal  trajectory.  Modal  control 
is  not  a  new  concept.  Wiesel  and  Shelton  (11)  proposed  a  modal  control  algorithm 
in  1983  to  control  a  satellite  in  an  orbit  about  the  exterior  Earth-Moon  Lagrange 
point  (Xi).  This  study  differs  from  theirs  in  that  we  will  use  impulsive  maneuvers 
to  control  the  satellite  instead  of  the  continuous  low-thrust  control  they  proposed. 

In  order  to  lead  up  to  the  problem.  Chapter  II  will  explain  the  dynamics  of  the 
restricted  three-body  problem,  which  is  a  Hamiltonian  system.  It  is  that  problem  on 
which  we  will  use  the  theory  of  Floquet.  The  modes  will  be  found  and  explained,  and 
then  a  controller  will  be  devised  using  impulsive  thrusters.  A  maneuver  is  performed 
when  the  unstable  mode  grows  beyond  a  certain  limit  —  the  "threshold”.  The 
direction  of  the  maneuver  is  determined  by  the  projection  of  the  imstable  mode  onto 
the  three  momentum  states.  We  will  run  simulations  to  verify  that  the  control  action 
worked  and  evaluate  the  control  energy  expended.  We  will  then  add  perturbations 
due  to  the  moon  and  the  eccentricity  of  Earth’s  orbit.  Simulations  will  be  run  to  see 
what  effects  these  perturbations  have  and  if  the  controller  still  works. 


II.  Background 

The  satellite  we  wish  to  control  is  not  in  orbit  about  the  earth,  nor  is  it  entirely 
in  an  orbit  about  the  sun.  It  is  in  what  is  called  a  ihrte-hody  orbit  where  the  satellite 
is  the  third  body  in  orbit  about  the  other  two  (called  primaries,  or  the  primary 
and  secondary  bodies).  The  three-body  problem  is  an  old  one,  having  been  studied 
by  Euler,  Lagrange  and  others  as  early  as  1772  (see  (7:  4)).  We  will  describe  the 
“restricted  three-body  problem”  and  find  equations  of  motion. 

2.1  Restricted  Three  Body  Problem 

Imagine  two  massive  bodies,  each  und«  the  gravitational  infiuence  of  the  other, 
rotating  about  their  common  center  of  mass  in  circular  orbits.  A  third  body,  assumed 
to  h  negligible  mass  —  it  does  not  affect  the  motion  of  the  primary  and  secondary 
bodies  —  orbits  under  the  infiuence  of  both.  Its  motion  is  affected  only  by  the  two 
primaries.  The  description  of  the  motion  of  this  third  body  is  the  solution  to  the 
restricted  three-body  problem  (7:  8). 

The  earth  and  sun  with  a  satellite  at  relatively  large  distance  firom  either  is  an 
example  of  a  system  for  which  the  restricted  three-body  problem  (RTBP  or  three- 
body  problem)  is  a  good  first  approximation.  The  sun  is  the  primary  (designated 
mi),  the  earth-moon  system  is  the  secondary,  mj,  and  the  units  are  chosen  so  that 
the  sum  of  their  masses,  the  distance  between  them,  and  the  gravitational  constant 
are  all  unity.  This  leads  to  the  fact  that  the  angular  frequency,  u,  of  their  motion 
(and  therefore  our  rotating  reference  frame)  is  unity  as  well  (see  figure  2.1). 

Defining  the  parameter  /i  =  mj,  we  have  mi  =  1  —  and  the  center  of  mass 
of  the  system  (remember  that  the  third  body  has  negligible  mass)  is  found  to  be 
a  distance  fi  firom  mi  leaving  1  —  p  between  center  and  mj.  We  set  our  origin  at 
the  center  of  mass,  with  the  positive  x  axis  pointing  toward  mi.  The  z  axis  points 
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Figure  2.1.  The  Three  Dimensional  Restricted  Three-Body  Problem 


upward,  and  the  y  axis  completes  a  right-handed  coordinate  system.  The  finune 

A 

rotates  with  d  =  Ifc. 

A 

The  third  body,  m3,  the  reference  satellite,  has  position  vector,  f=  xt-|-yj+xe 
in  rectangular  coordinates.  Differentiating  with  respect  to  time,  we  obtain 

•  A 

f  =  ii-i-yj  +  ih-f-wxr 

=  (*-y)»  +  (»  +  *)i  +  **- 

The  kinetic  energy  per  unit  mass  is 


T  =  if.f 


(2.2) 


The  potential  energy  per  unit  mass  is 


wlieie 


Pi  =  J(* -/*)*  + If* +  »> 


p*  =  +  +  +  (2.5) 

But,  lemembeiiiig  that  G  is  equal  to  one  in  otur  chosen  units,  and  substituting  the 
parameter  p  for  the  masses, 


I  M 

Pi  Pi 


The  Lagrangian  for  this  system  is 
L  =  T-V 


=  5  [(*-»)*  +  (»  +  *)*  +  **]- +  £  • 
^  i  Pi  Pi. 


Now  from  the  definition  of  generalised  momenta,  pi,  =  st',  we  have 

Oqh 


Pm  = 

ai  . 

(2.9) 

Py  = 

d£ 

(2.10) 

J>.  = 

If 

(2.11) 

and  we  see  by  comparing  equations  2.9-2.11  to  equation  2.1  that  the  generalised 
momenta  are  the  components  of  the  inertial  velocity  vector. 

n 

The  definition  of  the  Hamiltonian,  ^  Ph9k  —  ^  eliminating  qs's  in 
favor  of  pk'a,  yields: 

^  =  5(pi+pJ+i»J)+p-y-iv»- 

*  [  Pi  Pi. 
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HaiiultaiB*i  equAltoas  of  motioa  an  gmn  bj 


8H  .  8H 
9>  =  -sr>  Pfc- 


(J.1S) 


■o  perfanning  the  neceisaxy  partial  diflEierentiatkm  leada  to  the  equations  of  motion 


X 

y 

X 

1^ 

fv 

Pm 


= 

=  JV-* 

=  Pm 

=  IV  + 

=  -Pb  + 


A  A 

(1  -  <«)»  .  w 
(4  rf. 


(i-»i)»  ,  iti 

A  A 


(a.i4) 

(a.16) 

(J.16) 

(2.17) 

(2.18) 

(2.12) 


Now  let  the  state  vector  for  the  system  be  defined  by 


-(»»y,», Pm, Pg, Pm) 


(2.20) 


then 


=  (i,y,i,Pm,Pf,PM)- 


Introduce  the  matrix  Z,  which  is  defined  by 


■(:.:) 


(2.21) 


(2.22) 


where  /s  is  a  3  x  3  identity  matrix  and  0  represents  a  3  x  3  null  matrix.  Then  we 
can  write 


(2.23) 
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which  is  a  lestatement  of  Hamilton*!  equations  of  motion  in  state  vector  form.  This 
formulation  is  general  in  that  it  can  be  used  for  any  Hamiltonian  system. 

fB.2  Periodic  Orbits 

The  third  body  in  our  problem  is  in  orUt  about  the  primaries.  We  would 
like  to  find  out  what  its  orbit  looks  like  in  order  to  predict  and  control  the  state 
of  the  satellite.  If  this  were  the  two-body  problem,  the  answa  would  be  quite 
simple.  AH  two-body  orbits  are  conic  sections:  dOipses  for  relatively  low  energy 
orbits,  hyperbolas  for  the  high  energy  orbits,  with  parabolic  orbits  as  the  limiting 
case  in  between.  All  elliptic  two-body  orbits  axe  planar  and  periodic.  Given  a  period 
of  r,  the  state  vector  at  time  r  -I-  <o  i>  identical  to  the  state  vector  at  time  to> 

Not  all  the  orbits  in  the  three-body  problem  are  periodic,  but  periodic  orbits 
do  exist.  They  are  interesting  because  we  can  make  predictions  about  the  behavior 
of  periodic  orbits  over  long  periods  of  time  that  cannot  be  made  about  non-periodic 
orbits.  Finding  a  periodic  orbit  is  a  matter  of  guessing  initial  conditions  to  the 
state  vector  and  then  integrating  the  equations  of  motion  through  one  period.  The 
difiPerence  between  the  initial  state  vector  and  the  final  state  vector  is  the  error, 
e  =  z(r)— e(0).  The  initial  conditions  (and  the  period  if  necessary)  axe  then  adjusted 
until  the  error  vanishes.  In  order  to  make  corrections  to  the  initial  conditions,  the 
matrix,  $(r,  0),  needs  to  be  found.  $  is  the  differential  state  transition  matrix,  which 
relates  small  changes  in  the  state  according  to 

Sx{t)  =  m,to)Xito)  (2.24) 

where  Sx  =  X(t)  —  X(to).  This  is  a  linear  equation.  It  is  only  valid  in  a  region  of 
the  phase  space  where  X(t)  is  relatively  dose  to  X{to),  as  it  will  be  if  t  =  r  and 
to  =  0  and  provided  we  have  chosen  a  set  of  initial  conditions  near  a  periodic  orbit. 
$  is  found  from  the  variational  equations  of  the  system,  introduced  in  section  2.4. 
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t.S  The  Halo  OrM 


To  find  the  <»bit  fox  thii  itndy  we  fint  fonnd  the  value  for  the  dimenaioaleM 
parameter,  /i.  Let  mi  =  m0,  the  mam  of  the  inn,  and  mj  =  +  m0,  the  ram  of 

the  masaef  of  Earth  and  the  moon,  respectively.  With  a  little  algebra,  it  can  be  seen 
that 

„  =  _!!!_ 

•ni  +m> 

SO  we  obtain 

H  =  — —  =  3.040367143  x  lO"*.  (2.26) 

+  ffi^  +  m® 

Next,  we  looked  at  the  orbit  firom  ISEE-3.  Richardson  (6)  gave  approximate 
values  for  that  orbit.  We  transformed  the  approximate  initial  conditions  into  our 
own  reference  frame  (see  figure  2.1)  and  were  able  to  converge  to  a  periodic  orbit  Iqr 
iteratively  correcting  those  initial  conditions.  The  initial  conditions  we  converged  on 
were 

'  -.9916251461964399 

0.0 

-.0006706478525000000  ,  , 

X{U)  =  <  (2.27) 

0.0 

-.9818296716854701 

0.0 

\  / 

These  are  neither  exact  nor  unique.  Any  point  on  (or  very  near)  the  periodic  orbit 
could  be  used  as  initial  conditions,  and  will  yield  the  same  trajectory.  The  period 
was  found  to  be 

r  =  3.0596432056926  TU  (2.28) 

given  in  canonical  time  units  (described  below),  or  about  177.86  days.  The  first  three 
elements  of  the  state  vector  (x,  y  and  z)  are  in  distance  units  and  the  last  three  (pb,  Py 
and  p,)  are  in  distance  units  per  time  unit.  A  distance  unit  in  this  problem  is  an 
Astronomical  Unit  (AU)  which  is  approximately  1.495978  x  ^  1*  kilometers.  A  time 
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unit  if  —  timet  the  period  of  the  rotating  reCsrenoe  frame.  Since  that  period  it  one 
2ir 

year  for  thit  tyttem,  one  time  unit  it  about  2  sumtht  (58.132356144  dayt).  'Diete 
are  referred  to  at  canonical  unitt.  Referring  to  figure  2.1  for  the  reference  frame,  we 
tee  that  the  initial  oonditiont  place  the  tatdHte  about  one  hundredth  of  the  way  to 
the  tun  from  Earth,  tlightly  below  the  ecliptic  plane.  Since  the  orbit  it  not  planar, 
vituaJitation  can  be  a  rather  difilcnlt  problem.  Figure  2.2  thowt  the  projectiont  of 
the  orbit  onto  the  principle  (X-Y,  Y-Z,  and  X-Z)  planet  and  a  ‘^hree-dimentkmal* 
view  at  well. 

Thit  orbit  it  hereafter  referred  to  at  the  periodic  orbit  and  it  the  reference 
orbit  uted  in  thit  ttudy. 

£.4  Floquei  Theory 

A  periodic  orbit  it  a  tpedal  tdlution  to  the  equationt  of  motion  for  the  tyttem, 
jutt  at  an  equilibrium  point  it  a  tolution.  There  are  many  problemt  that  have  been 
tolved  by  lineariting  equationt  of  motion  about  an  eqnilibrinm  point  to  obtain  a 
conitant  coefficient  linear  tyttem,  and  in  much  the  tame  way,  we  Bneariae  about 
a  periodic  orbit  to  obtain  a  time  periodic  linear  tyttem.  A  Hamiltonian  tyttem  it 
linearized  by  taking  the  partial  derivative  of  the  equation!  of  motion  with  retpect  to 
the  state  vector. 

ex  ^  S>H 

dX  ~  ^dX* 

We  d^ne  A{t)  at  the  matrix  of  second  partial  derivatives. 

fntf 

^(‘)  =  (2-30) 
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Fignie  2.2.  Four  RepietentatioiiB  of  Halo  Orbit 


The  left  nde  of  Eq.  2.29  then  beomnet 


dm  dm 
dm  djt 

?i  ?t 

dm  dy 

?i 

dm  dy 

d^  d^ 

dm  dy 

^  % 

dm  dy 

dp,  dp, 
dm  dy 


dm  d^ 
d*  dp, 

dm  dp, 

?i  M. 

dx  dp, 

d^  dn 

dx  dp, 

?h.  ?h. 

dx  dp, 

dp,  dp, 
dx  dp. 


dm  dm 
dpt  dp, 

^  ?L 

dp^  dp, 

^  ?!L 

dpt  dp, 

d^  d^ 

dpg  dp, 

?h.  ?ii 

dpg  dp, 

dp,  dp, 
dpg  dp. 


(2.31) 


and  when  the  partial  derivatiYea  are  taken  for  the  BTBP  in  question,  A(t)  is  given 
by 

f  0  1  0  1  0  0  1 


-1  0 


A(t)  = 


0  1  0 
0  0  1 


A41  Ao  A4S  0  10 

Ail  Afi  Ais  —10  0 
Ail  Aii  Am  0  0  0 


(2.32) 
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•^1 

Att 

Au 


1-M  ,  3(l-/t)(g-M)*  M  . 

/»? 

>l,j  =  3»(»  -  m)(1  -  m)  _  SyM*  -  M  +  1) 
^  3*(* -/»)(!- m)  3*m(*-M  +  1) 

l-/>  3(l-;ty  H  ifiy^ 

A  ^  ri  A 


=  A«  = 


3|fx(l  -  m)  3|f*/» 


l-/>  .  3(1 M  .  3/tz* 

/»!  A 


(2.33) 


The  Imearuation  of  the  equations  of  motion  (2.23)  takes  the  form 


6x  =  A{t)Sx 


(2.34) 


Eq.  2.34  is  called  the  variational  equations  for  the  system.  The  state  transition 
matrix,  #(t,  to),  also  obeys 


*(»,«.)  =  A(0*(t.t.)  (2.36) 

By  integrating  Eqs.  2.23  and  2.35  (the  equations  of  motion  and  the  equations  of 
variation)  over  one  period  we  obtain  ^(r  +  to,  to),  the  monodromy  matrix.  This  is 
the  i  matrix  used  to  find  the  orbit  in  section  2.3.  #  has  no  obviously  meaningful 
structure,  and  it  will  not  be  given  here,  but  its  decomposition  into  eigenvalues  and 
eigenvectors  is  given  in  section  2.6. 
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We  leaxn  Cram  Fbqaet  (9)  that  the  lolutioa  to  the  vaaatkmal  equatioiu  can 
take  the  fonn  of  an  eigenirahie-eigenvector  proUem 

[#(r  +  to.to)-Ai/]/i  =  0  (2.36) 

where  Xi  axe  the  eigenwdues  and  ft  axe  the  corxefp<niding  ei^nvecton  of  the  mon- 
odxomy  matrix.  The  matrix  of  eigenveckm,  F,  will  be  important  becanie  it  proridef 
the  transfermatiQn  to  our  modal  ifaxiablei  (see  Chapter  m).  The  Poincatd  exponents, 
axe  defined  Iqr  Mj  3=  ^In  Aj,  so  the  eigenvalues  for  4^(r  +  to>to)  can  be  expressed 
as  Aj  =  e*^.  The  sdution  to  this  problem  has  txaditioaaUy  been  written  as 

m,  to)  =  (2.37) 

where  J  is  the  (usually)  diagonal  matrix  of  Poincard  e]q>onents,  but  this,  of  course, 
can  only  be  done  if  the  time  periodic  matrix,  F  =  ...)  is  invertible.  Assuming 

for  now  that  it  is  invertible,  we  can  differentiate  Eq.  2.37  to  get 

*{t,  to)  =  Fity^*-*^^F-\to)  +  F(t)  Je^<*-‘»)F-'(fo)  (2.38) 

and  substituting  the  right  hand  side  of  Eq.  2.35  into  Eq.  2.38  and  rearranging  we 

get 


FeJ(*-*»)F-\to)  =  A{t)F{t)e^^*-*^^F-\io)  -  F(t)Je-^<*-‘*>F-'(fo)  (2.39) 

Noticing  a  common  factor  in  Eq.  2.39  and  dropping  the  notation  for  functional 
dependence  on  time,  it  simplifies  to 

F  =  AF-FJ  (2.40) 

Equation  2.40  can  be  used  to  propagate  the  F  matrix  through  one  period. 
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t.S  EeUnded  Ei§enveetora 

There  are  cases  where  F  is  singolar,  aad  thereCrae  has  no  inverse.  In  &ct,  that 
is  the  case  with  the  halo  orbit  we  have  chosen  to  use.  The  #  matrix  has  a  repeated 
Poincar6  exponent  of  sero  (repeated  dgmivalue  of  nnity),  with  a  dependent  eigen¬ 
vector.  Since  F  =  any  eigenvector,  fi,  which  is  not  linearly  independent 

of  the  others  causes  F  to  be  rank  deficient  (singular).  In  our  case,  the  F  matrix  had 
rank  5  instead  of  6. 

The  generalised  eigenvector  presented  a  solution  to  our  problem.  For  the  setup 
of  the  generalised  eigenvector  problem,  let  ns  return  to  the  predecessor  of  £q.  2.37, 
which  has  no  F~^  term  (no  requirement  on  F  to  be  invertible). 

♦(t,  io)F(to)  -  F(t)e^t*-**J  =  0  (2.41) 


Note  that  the  takes  the  place  of  A,/  from  Eq.  2.36.  Let  us  use  the  Jordan 

form  of  the  2x2  block  from  the  J  matrix  which  contains  the  repeated  sero  value  as 
an  example  (the  Jordan  form  dictates  the  1  on  the  superdiagonal); 


where  we  have  let  t  =  r  and  to  =  0.  We  now  write  it  in  column  form  and  drop  the 
time  indices  for  convenience. 


'  ’ 

1  r 

0  1 

/ 


(2.43) 


Now  we  multiply  the  terms  above  by  columns,  which  leads  to 


-  {/xI(t/i  +  /i)}  =  0  (2.44) 
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SSnce  theie  two  m^rioet  axe  equivalent,  tkeir  cohnnm  ue  eqaivalent  ud  we  can 
write 

{♦  -  /)/.  =  0  (2.45) 

and 

{♦  -  /}/.  =  r/.  (X45) 

Equation  2.45  is  jnat  a  leatatement  of  Eq.  2.36  where  the  Poincai4  exponents  axe 
aero  =  0)  so  =  1.  In  order  to  generaliae  Eq.  2.46  we  write  it  as 


(♦  - 

where  it  the  repeated  eigenvector  and  ia  the  new,  generahied,  extended 
eigenvector.  Note  that  the  extended  eigenvector,  in  order  to  satisfy  Eq.  2.47,  can  be 
scaled  by 

/L  =  /-  +  “fee,  (2  «) 

and  it  is  therefore  as  non-uniqne  as  a  normal  eigenvector  (which  can  be  scaled 
/Lw  =  <x/oM  —  see  (10)).  We  will  now  use  Eqs.  2.36  and  2.47  in  continuing  onx 
example  to  calculate  the  eigenvectors  for  the  #  matrix.  First,  the  normal  eigenvalne 
equation  (Eq.  2.36)  is  used,  jpving  us 


fi 


(2.4#) 


The  vector  is  repeated  as  expected.  A  generaliaed  eigenvector  is  constructed  using 
Eq.  2.47  giving  us 

(x 

ftmt  =  span  I 

where  a  is  arbitrary  and  can  be  chosen  to  be  zero  (a  must  be  aero  if  the  new  vector 
is  to  be  perpendicular  to  the  old  one). 


(2.50) 
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2.6  ReavUa  of  tke  Bigenmilue-Eigenvecior  Prohlem 

The  eigenvalues  of  the  monodzomy  matrix,  $(r,0),  for  the  halo  orbit  were 
found  and  converted  into  Poincare  exponents.  The  Poincare  exponents  are  as  follows 

-2.4373955014588 
2.4373955019038 
-.026092033595943t 

Ui  =  < 

.026092033595943> 

0.0 
0.0 

\ 

The  first  two  exponents  are  real,  the  next  two  imaginary  and  the  last  ones  will  be 
referred  to  as  the  zero  modes.  The  software  which  produced  these  exponents  (the 
eigenvalue-eigenvector  package  in  the  IMSL  library  called  “devcrg”)  did  not  give 
them  to  us  in  this  order,  so  we  regrouped  them  for  ease  of  discussion.  The  order  is 
quite  arbitrary  as  long  as  each  exponent  is  associated  with  the  proper  eigenvector. 
It  should  also  be  noted  that  the  software  did  not  give  exactly  zero  for  the  zero 
exponents.  As  Wiesel  and  Pohlen  stated  (10:  4),  "It  is  up  the  worker  in  this  area  to 
recognize  this  case,  and  substitute  the  correct  values.” 

The  Poincare  exponents  are  constant  characteristics  of  the  orbit  in  the  three- 
body  problem.  The  negative  real  exponent  represents  the  stable  mode.  The  positive 
real  is  unstable.  Both  of  the  imaginary  exponents  represent  oscillatory  modes.  The 
zero  exponents  represent  only  a  different  epoch  of  the  same  orbit,  and  are  neither 
stable  nor  unstable.  Within  the  linear  region,  plots  of  the  first  two  modes  against 
each  other  look  like  a  saddle  point.  Plots  of  the  second  pair  are  circles  and  the  last 
two  remain  constant  in  time.  This  analysis  of  the  modes  is  the  basis  of  the  modal 
control  schemes  presented  in  Chapters  IV  and  V. 

Associated  with  each  of  the  exponents  is  an  eigenvector,  the  combination  of 
which  is  the  modal  matrix,  F{t),  which  was  used  in  the  similarity  transformation 


described  in  section  3.1.  The  mntriz  itsdf  is  Toy  cambets(»ne  and  it  is  pesiodic  (ior 
the  periodic  orlnt),  which  means  that  it  is  different  at  every  time,  t,  bnt  at  the  initial 
time  it  is  given  by 


m = { 


-2.374442084986946JS  -  01 
2.041548870256156£  -  01 
8.904551229102137.^  -  03 
-9.274729220416125E  -  01 
1.938588253107533£  -  01 
6.383160385046611£  -  02 


-2.374442084873532^  -  01 
-2.041548869917724^  -  01 
8.904551230670204£  -  03 
9.274729220687803£  -  01 
1.938588252272955£  -  01 
-6.383160385939196E  -  02 


-1.340603252397338^  -  02 
4.197004838119940£  -  09 
4.606002014311875.E  -  01 
-3.420895405951368£  -  09 
2.049158967689975E  -  01 
-5.203976092078526^  - 10 


-1.407327854914285£  -  27 
6.549857854528899£  -  01 
-2.845632307073551^  -  09 
-5.40614954061U)01£;  -  01 
-7.710108270570280J5  - 10 
-1.562251589378517B  -  01 


o.ooooooooooooooof; + oo 

7.415056092412337£  -  01 
O.OOOOOOOOOOOOOOOE  +  00 
-5.936711168454215£?  -  01 
O.OOOOOOOOOOOOOOO.S  +  00 
3.126084395650520£  -  01 


-2.217844788884213£;  +  01 
2.871181813841697£  -  08 
-1.495894728274239^  +  02 
-2.530885149026978£:  -  08 
5.398535426318291£  +  01 
-1.161680356454845^  -  07 


The  last  column  is  the  ertended  eigenvector. 


2.7  Summary 

We  have  equations  of  motion  and  of  variation  for  a  particle  in  a  three-body 
orbit.  Using  ISEE-3  data,  we  have  found  a  periodic  orbit  to  put  our  satellite  into. 
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We  kaTe  eeidoaled  the  Poincai^  ezpoMttte,  ead  naiiig  extended  e^eatTeeton  «e  hme 
found  an  initial  eigenvect<»  matzix,  F(0),  which  has  foil  nmk. 


III.  Theory 


S.l  Canameal  l^mafortnations 

In  order  to  do  perturbation  work  on  Hamiltoninn  lyitenu,  it  is  often  quite 
desirable  to  transform  the  system  into  another  set  of  variables  which  simplify  the 
construction.  If  the  transformation  preserves  the  Hamiltonian  structure  of  the  sys¬ 
tem,  it  is  called  a  canonical  transformation.  It  is  sometimes  possible,  in  fact,  to 
transform  to  a  set  of  variables  in  which  the  Hamiltonian  of  the  system  vanishes  iden¬ 
tically.  Such  a  transformation  solves  the  problem,  and  the  variables  are  ideal  for 
use  in  perturbation  work  since  the  solved  part  of  the  problem  no  longer  encumbers 
the  statement  of  the  dynamics.  We  were  not  so  ambitious  as  to  try  and  solve  the 
problem  because  the  RTBP  cannot  be  solved  in  closed  form  —  see  Saebehely  (7:  7). 
Our  goal  was  to  transform  the  problem  into  a  set  of  variables  that  allowed  the  funda¬ 
mental  modes  to  be  independent  of  one  another.  If  $  were  diagonal,  the  subspaces, 
or  modes,  of  #  would  be  independent,  and  that  is  what  led  us  to  the  eigenvalue- 
eigenvector  problem  in  Chapter  2.  This  set  of  variables  is  called  modal  variables,  and 
will  be  designated  q,-.  Wiesel  and  Pohlen  (10)  showed  how  to  do  canonical  transfor¬ 
mations  in  Hamiltonian  systems  for  the  non-degenerate  case  and  several  degenerate 
cases  which  appear  in  the  RTBP.  The  transformation  to  modal  variables  is  done  by 
noting  that  F  is  the  modal  matrix,  which  transforms  the  i  matrix  to  diagonal  form 
according  to 

A  =  (3.1) 

This  is  nearly  a  restatement  of  Eq.  2.41  where  e*'^,  the  exponential  of  the  Jordan 
form,  takes  the  place  of  A,  the  diagonal  matrix.  We  define 

Sx(t)  =  X{t)  -  (3.2) 
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to  be  a  imall  Taiiation  in  the  state  vector  at  a  partknlar  time.  A^(t)  is  the  state 
vector  for  the  periodic  orbit,  and  Jf(t)  is  the  state  vector  for  a  nearby  orbit,  not 
necessarily  periodic.  A  nearby  orbit  is  one  whose  state  vector  is  dose  to  that  ci  the 
periodic  orbit.  Note  that  Sx{t)  is  a  vector  of  length  six  called  the  variation  of  X. 
Then 

<*(«)  =  (3.3) 

where  F{t)  is  the  modal  matrix  at  that  particular  time  so  the  transformation  is 
simply  a  matter  of  solving  for  ii{t) 

Mi)  =  F-\i)SMt)  (3.4) 

The  modal  variables  approach  aero  as  the  state  vector  approaches  the  periodic  orbit. 
The  value  of  each  mode  changes  with  time,  ipving  a  plot  of  the  modes  a  characteristic 
shape,  and  a  characteristic  e*folding  time,  according  to  its  Poincare  exponent.  The 
e-folding  time  is  the  amount  of  time  required  for  the  mode  to  change  by  a  factor  of 
e.  Plots  of  the  mode  shapes  will  be  shown  in  the  next  chater. 

3.2  Fourier  Repreaentationa 

In  order  to  control  the  satellite  with  modal  control,  it  was  necessary  to  evaluate 
the  modal  variables  at  arbitrary  times  during  the  simulation.  Since  F(t)  varies  with 
time  (it  is  periodic),  but  not  according  to  a  simple  function,  there  had  to  be  a  way 
to  obtain  values  in  double  precision  for  the  6x6  matrix  at  arbitrary  times.  This 
was  accomplished  by  using  a  truncated  Fourier  series.  Each  member  of  the  periodic 
matrix  F  was  stored  as  a  table  of  sine  and  cosine  coefficients  so  that  any  given  value 
of  the  matrix,  /»,  was  given  by 

/«  =  S  co8(j(ot)  +  ^*k  an(ka;t).  (3.5) 

iso  k=i 
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This  sum  could  be  evaluated  at  any  time  during  the  simulation,  which  allowed  access 
to  any  F(t)  without  re-integrating  Eq.  2.40.  The  values  for  the  periodic  orbit 
were  stored  and  retrieved  in  the  same  way  alter  they  >vere  obtained  by  int^rating 
the  equations  of  motion  and  variation  in  the  beginning  of  the  project.  Fm  all  the 
calculations  done  in  this  study,  25  coefficients  were  used  and  each  coefficient  was 
calculated  and  stored  to  13  decimal  places.  This  provided  all  the  accuracy  needed 
for  the  calculations,  but  there  were  instances  where  the  truncation  of  the  Fourier 
series  represmtations  became  apparent.  Chapter  IV  shows  some  graphs  (see  figure 
4.6)  that  may  indicate  anomalies  due  to  this  approximation. 
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IV.  Controllers 


4.1  Tke  Halo  OrM 

The  hak>  orbit  is  inherently  nnstnUe.  Using  initial  conditions  whidi  wetc  cor¬ 
rect  to  within  ten  centimeters  in  position  and  a  lihe  amount  in  momentum,  the 
satellite  left  the  nominal  orbit  completely  befine  it  had  completed  its  &nrth  resolu¬ 
tion  (see  Figure  4.1).  When  evaluating  the  relative  stnbiKty  of  the  orbit  like  this,  care 
must  be  taken  that  enough  integration  steps  ate  used.  If  the  steps  are  not  suffidentbr 
dose  together,  the  errors  in  the  integration  can  ^  ^  inherent  instabilities  in 

the  orbit.  The  integration  package  used  was  a  fourth  order  predictor-corrector  (9).  It 
required  about  2000  integration  steps  per  orbit  to  assure  that  the  errors  were  smaller 
than  the  instability  so  the  true  instability  was  revealed.  The  divergence  shown  here 
resulted  from  a  tiny  error  in  the  initial  conditions,  which  grew  (exponentially  in  the 
unstable  modal  direction)  until  the  satellite  literally  coasted  off  into  space. 

4.1  The  "CheaUr"  CotUrotter 

The  "Cheater”  Controller  was  devised  as  a  way  to  stabilise  the  orbit  through 
several  periods  so  that  the  long  term  behavior  of  the  orbit  and  its  modes  could  be 
analysed  and  visualised.  The  controller  worked  by  resetting  to  sero  the  unstable 
mode  any  time  that  it  reached  a  threshold  value.  This  was  accomplished  by  sub- 
racting  from  the  state  vector  the  fraction  of  the  corresponding  unstabk  eigenvector 
equal  to  the  magnitude  of  the  threshold. 

~  ~~  */(^)i»iJtaM«  (^‘l) 

where  e  is  the  threshold  value.  This  simple  concept,  while  not  a  realistic  contrdler, 
validated  the  concept  of  modal  analysis  by  demonstrating  that  the  modes  were  in¬ 
dependent  of  one  another  since  a  change  in  the  state  vector  in  the  direction  of  the 
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Fignxe  4.1.  X  —  Y  Projection  of  Orbit  Showing  that  the  Satellite  does  not  Remain 
Close  to  the  Nominal  Orbit.  The  Earth  is  at  About  X  =  -1,  Y  =  0. 
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mutable  mode  did  not  affect  tbe  other  modes.  Using  this  oontioller  «e  vere  aUe  to 
propagate  tbe  mint  Cor  as  long  as  desired,  ence  it  was  not  allowed  to  diverge,  'nis 
was  important  because  it  turned  out  that  tbe  oscillatory  mode  bad  a  period  of  about 
seventy  times  tbe  orbital  period.  Figure  4.2  shows  tbe  two  imaginary  (oscillatory) 
modes  plotted  against  each  other  for  about  99  orbital  periods  (99r).  During  a  sim¬ 
ulation  of  only  a  few  times  r,  tbe  oscillation  completes  only  a  portion  d  tbe  circle 
seen  in  tbe  figure.  The  initial  conditions  fer  this  simulation  were  perturbed  in  tbe 
direction  of  one  of  tbe  oscillatory  modes  by  a  fector  of  10~^  so  that  tbe  mode  did 
not  just  hover  about  sero.  Tbe  orbit  was  controlled  using  tbe  ‘‘cheater*  controller 
with  the  threshold  set  to  10~^. 

The  reason  this  controller  is  unrealistic  is  that  tbe  entire  state  vector  was 
repositioned  during  tbe  ‘4nanenver*.  In  reality,  it  is  only  possible  to  change  tbe 
momentum  states  using  thrust.  The  concept  of  changing  a  position  state  instantap 
neously  corresponds  to  teleportation.  While  perhi^u  desirable,  it  is  not  possiUe,  so 
we  were  then  forced  to  look  at  a  more  feasible  alternative. 

4.S  Maneuver  Optimization 

A  realistic  controller  is  an  order  of  magnitude  more  complex.  Tbe  cheater  had 
the  advantage  of  moving  tbe  state  vector  in  the  direction  of  an  eigenvector,  affecting 
only  one  mode.  Whatever  we  do  with  an  actual  maneuver  will  affect  tbe  state  vector 
in  a  somewhat  more  unpredictable  way  than  our  cheater  did.  If  we  are  to  obtain 
the  desired  response  in  one  mode,  we  will  have  to  sacrifice  some  response  in  the 
others.  This  amounts  to  the  statement  of  an  optimisation  problem,  and  indeed,  in 
this  chapter  we  attempted  to  optimise  control  by  two  separate  strategies.  The  first 
is  described  in  this  section,  and  the  second  in  the  next. 
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Mode  Four 


Fignze  4.2.  Osdllstoiy  Modes  Plotted  Against  One  Another  for  99  Periods 
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of  tlie  state  sector,  SX^ 


Then 


=  {0  0  0  4).  %  (4J) 

since  the  momentum  states  can  be  aAscted  bj  an  impnkhre  maaensec.  Recall 
the  transformation  to  modal  saziables  from  Eq.  3.4. 

ij  =  F-'SX  (4.3) 

We  can  then  write  the  change  in  modal  sariaUes  as 

Sti  =  F'HX  (4.4) 


where  the  S  notation  is  used  here  to  denote  a  manenser  (a  change  in  state).  When 
we  write  each  of  the  cmnponents,  it  looks  like  this 


Where  F^  are  the  components  of  F~^  (and  not  the  inverse  of  the  components 
of  F).  Note  that  we  will  be  required  to  compute  the  actual  inverse  of  the  modal 
transformation  matrix,  F,  for  the  first  time.  Before,  we  had  always  been  able  to 
solve  the  linear  algebra  problem  for  the  modal  variables  by  Gaussian  elimination 
and  back  substitution  withcmt  computing  the  inverse  (see  Chiq>ter  m).  Since  Sifi  is 
the  variation  in  the  nnstabfo  mode,  we  axe  interested  in  the  elements  of  F~^  which 
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boUl  oontrilrate  to  the  vihie  of  Sih  aad  mnltiitly  the  pennkaUe  derneuU  of  SX. 
Theie  axe  jP{j4 ,  aad  This  constitutes  a  sia^  sodax  eqnatkm 

Stfi  =  FiJSpg  +  FiJSpn  +  F{JSpm  (4.6) 

Just  as  in  the  "cheatex”  oontxdlex,  we  would  hhe  to  the  unstable 

mode.  Mathematically, 

^VlAuit*d  —  ~~Vl- 

The  desixed  change  (maneuvex)  in  the  fixst  modal  vaiiable  is  to  negate  its  pxesent 
value.  This  constitutes  a  second  scalax  equation.  We  axe  now  leady  to  foxinally  state 
the  optimisation  pxoblem.  We  wish  to  minimise  the  fnd  used  fox  station-keeping. 
We  will  use  AV  pex  unit  time  as  the  fignxe  of  mexit  to  evaluate  the  contidlexs.  The 
contxollex  will  minimise  P  =  (AV)^  for  eacl  i  maneuver.  Stated  mathematically,  we 
will  minimise  Eq.  4.8  subject  to  Eq.  4.9. 

(4«) 

Svi+Vi=0  (4.9) 

Substituting  from  Eq.  4.6  fox  the  value  of  and  introducing  a  Lagrange  multiplier, 
A,  allows  us  to  wxite  the  pxoblem  in  one  equation. 

P  =  6pl  +  Sp^  +  Sfijl  +  A(Fi~4ip,  +  F{j^6p^  -I-  +  qi)  (4.10) 
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IkUag  putki  dflrimthras  of  P  with  icspeel  to  the  pemiiMUe  chugM  in  ike  state 
aad  setting  them  eqval  to  aero  yidds 


^  =  o  =  it,.+xFri  =>  ^  = 
j^=o=iffc+xjr;  =»• 
^=o=Mp.+Afs  =►  fi..=-jAjir; 


(4.11) 


and  with  respect  to  the  Lagrange  multiplier 


dP 

8X 


=  0  —  6rfi  +11%  Sifi  =  —ifi. 


(4.12) 


Note  that  Eq.  4.12  is  simply  a  restatement  of  the  constraint.  Ftom  the  right  hand 
side  of  Eq.  4.11  we  know  the  direction  of  our  maneuver  since  we  know  the  proportion 
in  each  of  the  momentum  directions,  but  its  magnitude  is  scaled  by  A,  which  we  do 
not  know.  In  order  to  solve  for  it,  let  us  rewrite  Eq.  4.11. 


Sp  = 


1. 

(  f-A 

^X,4 

Spy 

-"2^ 

Then  note  that  Eq.  4.6  can  be  rewritten  as 


(4.13) 


%  =  (n:i  fri  Friw 


Substituting  from  Eq.  4.13  for  the  value  of  Sp  we  have 

S’h  =  -iAK^M')'  +  )’  +  (^M  « 


(4.14) 


(4.1S) 
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We  know  from  tke  ooutnint  that  ^  s  91,  lo  lolTing  for  A 


A  = 


391 _ 

.  w;)* + + (ffi)*. 


(4.W) 


The  contnd  Kheme,  then,  conncted  of  evaluating  the  unitaMe  mode  until  a 
threshold  was  rmwhed.  Since  the  unstable  mode  grows  exponentially  with  time,  we 
knew  it  would  not  decrease.  A  realistic  impulsive  maneuver  was  planned 
to  the  mathematics  presented  in  this  section  which  optimised  the  direction  and  the 
magnitude  of  the  maneuver  to  sero  out  the  unstable  mode.  The  equations  of  motion 
were  then  prop<^ted  from  that  point,  with  initial  conditions  given  by  the  conditions 
at  maneuver,  plus  the  maneuver  (the  change  in  momentum  states),  until  the  unstable 
mode  again  reached  the  threshold. 


Jl^.\  ResuUa  of  the  ReaHaUc  Controller 

The  results  using  this  controller  are  displayed  in  graphical  form  in  the  following 
figures.  We  begin  with  a  medium  thxeshdd,  allowing  the  unstable  mode  to  reach  a 
value  of  10~^.  See  Figure  4.3.  This  threshold  was  chosen  because  its  results  axe  a 
good  compromise  between  low  cost  and  frequency  of  maneuvers.  A  lower  cost  can  be 
obtained,  but  it  requires  more  frequent  maneuvering.  The  average  normalised  cost 
of  the  maneuvers  approaches  a  level  value  after  about  40  time  units  which  is  a  little 
over  13  orbits  and  about  six  and  a  half  years.  The  simulations  in  this  section  were 
run  for  a  total  of  100  time  units-some  30  orbits  over  about  15  years.  The  average 
cost  for  maintaining  the  orbit  for  this  first  threshold  is  about  0.6  cm/sec  per  year. 
Each  of  the  maneuvers  is  about  0.4  cm/sec  in  magnitude  and  th^  occur  with  a 
frequency  of  just  less  than  once  per  orbit  (just  under  twice  a  year). 

Figure  4.4  shows  the  performance  of  the  same  controller  with  a  lower  thresh* 
old  of  10~*.  Here,  the  maneuvers  ate  happening  with  greater  frequency  (about 
double-average  is  1.7  maneuvers  per  orbit,  which  is  3.4  maneuvers  per  year),  but 
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Figure  4.3.  Average  Maneuver  Cost  as  a  Function  of  Time  (Canonical 
Threshold  =  10”^ 
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Figure  4.4.  Average  Maneuver  Cost  as  a  Function  of  Time  (Canonical  Units), 
Threshold  =  10~® 

the  magnitude  of  each  maneuver  is  much  smaller  (they  average  about  0.0037  cm/sec). 
The  total  average  cost  is  only  0.013  cm/sec  per  year-about  l/50th  what  it  was  to 
maintain  it  with  the  higher  threshold.  Now  let  us  look  at  the  behavior  of  the  modes 
during  this  control  process.  Figure  4.5  is  a  plot  of  the  two  real  modes,  correspond¬ 
ing  to  the  positive  real  Poincare  exponent  (unstable  mode)  and  the  negative  real 
exponent  (stable  mode).  The  general  shape  of  the  plot  may  be  a  bit  deceiving  at 
first.  The  diagonal  lines  are  drawn  where  the  maneuver  occurs.  The  unstable  mode 
is  brought  back  to  (approximately)  zero,  but  the  stable  mode  acquires  an  unknown 
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Figure  4.5. 


politive  or  negative  value  (remember  we  can  only  change  momentum  states,  so  the 
controller  affects  other  modes).  As  the  orbit  progresses,  the  stable  mode  approaches 
aero  while  the  unstable  mode  diverges.  This  creates  a  saddle  shi^  in  the  center 
of  the  plot,  as  shown  by  Figure  4.6.  The  spikes  in  the  figure  are  due  to  numerical 
errors  in  calculating  the  modes  (possibly  truncation  error  in  the  Fourier  coefficients). 
Larger  values  of  the  modes  are  less  affected  by  these  errors,  while  smaller  values  are 
nearly  obscured  by  them.  Figure  4.7  would  look  very  messy  if  all  the  dots  were 
connected  because  of  the  large  errors  near  the  origin,  but  as  a  series  of  dots,  it  still 
shows  the  same  progression  u.  the  real  modes.  The  unstable  mode  is  never  allowed 
to  exceed  10~*  in  absolute  value. 

It  would  be  interesting  to  see  what  would  happen  if  we  tried  to  reduce  the  num¬ 
ber  of  maneuvers  further  by  increasing  the  threshold.  We  could  allow  the  unstable 
mode  to  get  larger  before  imposing  control.  Figure  4.8  shows  that  the  cost  continues 
to  increase  as  time  progresses  for  a  threshold  of  10~*.  Although  the  first  maneuver 
happens  later  than  it  does  with  a  smaller  threshold,  they  come  more  frequently  after 
the  first  two  or  three.  This  is  not  hard  to  understand  if  you  look  at  the  plot  of  the 
real  modes  in  Figure  4.9,  which  shows  that  the  modes  are  really  no  longer  in  the 
linear  region  where  they  behave  like  a  saddle  p<^t  and  where  they  can  be  controlled. 
The  controller  still  does  not  allow  the  unstable  mode  to  exceed  the  threshold,  but 
this  is  insufficient  to  control  the  orbit  in  this  nonlinear  region.  One  interesting  thing 
about  this  control  scheme  was  that  the  orbit  itself  did  not  remain  the  same.  Figure 
4.10  shows  the  view  of  the  orbit  that  we  would  get  from  the  earth.  It  remained  a  halo 
orbit  for  about  three  and  a  half  periods,  and  then  began  to  change.  It  was  under 
some  control,  however,  because  it  failed  to  diverge  completely.  Its  apparent  direction 
reversed  as  the  plane  of  the  orbit  revolved.  The  worst  thing  about  this  situation  was 
that  the  satellite  crosses  the  sun’s  disk,  which  was  completely  unacceptable.  The  sun 
would  appear  at  the  center  of  the  plot,  about  1  inch  in  diameter.  The  cost  diverged 
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stable  Mode 


Unstable  Mode  ^  ^q*io 


Figure  4.7.  Plot  of  the  Two  Real  Modes  as  the  Realistic  Controller  Operates  Show¬ 
ing  the  Saddle  Behavior,  Threshold  =  10~* 
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Fignie  4.9.  Plot  the  Real  Model  for  am  Inraffident  Amount  of  Control,  Threahold 
=  10~*,  Showing  Nonlinear  Behavior  of  the  Model 
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Figure  4.10.  Plot  of  the  Y  -  Z  Projection  of  the  Orbit,  Threshold  =  10~*,  Showing 
that  the  Orbit  Changes  RadicaOj  Whik  Remaining  in  Vicinity  of  Lt 
(Origin) 
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Figure  4.11.  Performance  of  the  Realistic  Contrdler  as  a  FVmction  of  the  Con¬ 
troller’s  Thresh<dd  Value 

and  the  satellite  left  its  orbit.  We  concluded  that  the  unstable  mode  could  not  be 
allowed  to  reach  10~*  in  magnitude  if  we  wished  to  contnd  the  satellite. 

The  figure  of  merit  used  to  compare  the  contnd  schemes  is  the  average  ac¬ 
celeration.  Each  maneuver  is  a  velocity  change,  AV.  The  only  way  to  objectively 

evaluate  the  performance  of  the  controller  is  to  normalise  the  total  AV  by  the  time 
AV 

elapsed,  which  is  an  average  acceleration.  Figure  4.11  shows  the  perfnrmance 
of  this  controller  as  a  function  of  the  threshold  value.  We  can  see  that  less  maneu¬ 
vering  is  required  if  the  threshold  is  kept  lower.  The  station-keeping  cost  has  been 
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oolivcxted  to  naitt  Uiat  are  eaay  to  ui^entaad  (cm/aec  per  year)  for  tUs  plot,  bet 
tbe  boiiaontal  axis  is  left  in  tbe  <»i^mal  “iiKxial”  aaits,  wbidi  are  a  tzansformatioa 
(by  F~^)  of  canonical  nnits.  Each  point  on  the  plot  in  Figure  4.11  represents  a  mn 
(simulation)  of  the  system,  and  is  the  final  ralne  of,  for  instance,  Figues  4.3  and 
4.4. 


4-5  OpHmizaiion  of  Manenner  Epoch 


In  this  section,  a  slightly  different  centred  scheme  was  attempted.  Instead  of 
immediately  maneavering  whenever  the  unstable  mode  reached  a  threshold,  the  ideal 
time  for  a  maneuver  in  the  orbit  was  determined,  and  then  a  maneuver  was  executed 
at  that  time  regardless  of  the  instability.  The  optimal  time  for  a  maneuver  is  when 
the  projection  of  the  unstable  mode  onto  the  momentum  direction  is  maximum 
80  that  a  momentum  change  will  produce  mMrimnTn  change  in  the  unstable  mode. 
The  concept  of  optimising  the  maneuver  epoch  and  the  maneuver  direction  and 
magnitude  should  give  us  a  global  optimum.  It  remains  to  be  seen,  however,  whether 
that  global  optimum  is  suited  to  the  problem.  The  reason  it  might  not  be  is  that 
the  e-folding  time  of  the  unstable  mode  (time  required  for  it  to  grow  by  a  factor  of 
e)  is  seven  orders  of  magnitude  shorter  than  the  period  of  the  orbit.  This  results  in 
the  growth  of  the  unstable  mode  by  a  factor  of  e^  (approximately  =  1100)  during 
one  period. 


The  optimum  time  for  a  maneuver  is  when  the  maneuver  itsdf  would  be  small¬ 
est.  In  the  development  of  the  previous  controller,  we  minimised  the  objective  func¬ 
tion,  P(=  AV^),  for  whatever  time  in  the  orbit  that  the  threshold  happened  to  be 
exceeded.  Now  we  wish  to  find  a  global  minimum  for  AV  regardless  of  the  time  or 
the  value  of  the  unstable  mode.  In  order  to  find  this  global  minimum,  we  remove 
the  modal  variable  from  the  problem.  The  normalised  value  of  A  becomes  (see  Eq. 
4.16) 


A  = 


2 


(4.17) 
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■ad  ire  know  tkai  tlie  magaitade  of  tke  nl^awlaer  ie  grrea  bj 

AV  -  +  ^  +  «  (4.18) 

■o  eubatitutiag  foot  the  iralue  of  the  Sff»  from  Eq.  4.11 

AV  =  +  (-iwri)*  +  (4.W) 

When  we  inbetitnte  from  Eq.  4.17  for  tlm  valne  of  A,  and  perfbnn  the  obvioiu 
cancellation,  we  see  that 

AV  =  VA  (4.20) 

so  that  tnifiitnidiig  AV  u  equivalent  to  fniwmiMwg  the  nonnaUaed  A.  It  tuna  ont 
that  A  has  a  periodic  nature,  and  has  the  aame  period  aa  that  of  the  orbit.  Figure 
4.12  ahowa  the  value  of  the  normaUaed  A  aa  a  function  of  time  for  10  time  nnita 
(3.27  perioda).  Ita  value  flnctnatea  between  about  2.18  and  2.06,  with  a  minimum 
occurring  jut  once  per  orbit. 

Knowing  that  A  wu  periodic  and  that  it  had  a  minimum  during  eadi  orbit 
allowed  u  to  uae  that  fact  in  maneuvering  the  aatellite.  The  contnd  acheme  evaluated 
A  at  each  point  in  the  orbit,  found  where  it  had  ita  minimum,  and  then  executed 
a  maneuver  at  that  time.  The  maneuver  wu  computed  in  exactly  the  aame  way 
u  wu  done  in  aection  4.3.  The  performance  of  the  contrdkr  ia  ahown  gnq>hicany 
in  Figure  4.13.  Compariaon  of  Figure  4.13  with  Figure  4.3  ahowa  that  thia  method 
of  maneuvuing  ia  actually  a  little  leas  expenaive  (about  0.23  cm/aee  per  yeu)  than 
the  threshold  controller  with  the  threshold  at  10~^,  but  not  u  cheap  u  when  the 
threahold  is  moved  to  10~*.  Another  advantage  of  the  periodic  contrdler  here  is 
that  after  the  first  maneuver  (0.86  cm/sec),  all  the  othera  were  neariy  identical  in 
magnitude  (the  maneuver  magnitude  grew  slowly  from  about  0.07  to  0.17  cm/sec) 
and  direction,  and  they  were  all  performed  at  the  same  point  in  the  orbit.  See 
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Normalized  Value  of  Lagrange  Multiplier 
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Figure  4.13.  Average  Maneuver  Cost  as  a  Function  of  Time  Wliile  Maneuvering 
Once  per  Orbit  at  the  Optimal  Time 
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Figure  4.14.  X  -  Y  Projection  Plot  of  the  Maneuvering  Portion  of  the  Halo  Orbit 
Under  Periodic  Control 

Figure  4.14.  The  maneuvering  point  is  where  several  points  are  plotted  tc^ther  in 
the  middle  of  the  trajectory  shown.  The  vdodty  change  at  this  point  is  so  slight  as 
to  be  invisible  to  the  eye  even  at  this  scale,  which  shows  about  a  hundredth  of  the 
orbit. 

The  performance  of  this  controller  was  somewhat  disappointing.  While  the 
concept  of  the  optimum  maneuver  epoch  was  validated  (it  worked),  it  was  not  overall 
the  cheapest  way  to  control  the  satellite  in  its  orbit.  It  was  not  therefim  a  ^bal 
optimum  for  the  control  problem.  A  low  threshold  (say  10~*)  was  cheaper  by  an 
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order  of  magnitude.  If  A  had  had  a  greater  range,  so  that  there  was  a  point  in  the 
orbit  at  which  it  was  much  cheaper  to  maneuver,  then  this  optimum  would  have 
held  greater  promise  for  inexpensive  maneuvering.  As  it  is,  the  optimum  epoch  fi>r 
maneuvers  will  not  be  pursued  further. 

4.6  Behavior  of  the  Other  Modes 

While  it  is  not  strictly  on  the  critical  path  to  controlling  the  satellite,  the 
behavior  of  the  other  modes  under  control  is  interesting  to  note  and  will  therefore 
be  included  here  briefly.  It  was  noted  in  section  4.2  that  the  oscillatory  modes 
would  strictly  oscillate  about  the  origin  if  the  unstable  mode  could  be  independently 
controlled  within  the  linear  region  (see  Figure  4.2).  If,  however,  the  unstable  mode 
is  controlled  with  some  dependence  on  the  other  modes  (as  it  was  with  the  realistic 
controllers),  we  might  wish  to  know  what  the  the  crossover  effect  will  be. 

Let  us  look  at  Figure  4.15.  During  each  of  the  maneuvers  there  is  a  disconti¬ 
nuity  in  the  plot,  shown  with  a  straight  line  (more  or  less  horizontal  Hnes).  Between 
maneuvers,  when  the  satellite  is  drifting  freely,  the  imaginary  modes  tend  to  con¬ 
tinue  on  their  characteristic  circles  (nearly  vertical  arcs  on  the  graph).  The  period 
of  the  imaginary  mode  oscillation  is  so  long  that  only  a  short  segment  of  the  circle 
is  completed  before  the  next  maneuver.  The  value  of  the  imaginary  modes  appears 
somewhat  random,  but  it  increases  slowly  over  the  course  of  the  simulation. 

Figure  4.16  is  the  analog  of  the  previous  figure,  with  the  zero  modes  plotted 
instead  of  the  imaginary  modes.  When  the  satellite  is  drifting  freely,  the  value  of  the 
modes  remains  relatively  constant  (these  appear  as  vertices  on  the  graph),  and  when 
a  maneuver  is  performed,  their  value  changes  slightly  (the  whole  plot  is  on  the  order 
of  10~*)  in  the  direction  of  the  maneuver,  to  remain  there  until  next  maneuver. 
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Mode  Six 


ModeRve 

Figure  4.16.  Zero  Modes  Plotted  Together  While  Satellite  is  Under  Realistic 
Control 
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Swmnary 

It  has  been  demonstrated  that  the  satellite  in  the  three-body  halo  orbit  can 
successfdlly  be  controlled  (kept  near  the  periodic  orbit)  by  observation  and  manip¬ 
ulation  of  its  modal  states  with  little  attention  being  paid  to  its  actual  position  in 
space  or  velocity  at  any  given  time.  If  the  unstable  mode  is  carefully  watched  and 
never  allowed  to  exceed  a  pven  threshold,  the  satellite  will  remain  very  dose  to 
the  desired  periodic  orbit  as  long  as  the  threshdd  is  low  enough.  It  has  also  been 
shown  that  for  the  unperturbed  case,  this  form  of  control  is  very  inexpensive  in  terms 
of  total  maneuvering  fud  expended.  The  optimization  of  the  maneuver  itself  was 
superior  to  optimization  of  the  maneuver  epoch  in  total  station-keeping  cost.  The 
control  cost  could  be  made  very  low  by  sdection  of  a  smaller  threshold.  The  epoch 
optimization  controller,  though,  could  be  more  practical,  because  the  maneuvers  are 
nearly  constant  and  occur  at  regular  intervals  and  could  be  scheduled  and  calculated 
in  advance. 
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V.  Perturbationa 


The  saidlite  in  the  thiee-body  orbii  has  been  oontzdled  (see  Chiyiter  IV).  If 
Earth  were  in  a  drcolai  orbit  about  the  sun  with  no  moon  or  other  perturbaticms, 
our  contrcd  problem  would  be  s<dred.  It  is  the  purpose  of  this  chapter  to  introduce 
the  perturbations  to  the  system  due  to  the  motion  of  the  moon  about  the  Earth  (we 
will  call  this  effect  ‘^oon’*),  and  the  eccentricity  of  the  Earth’s  orbit  about  the  sun 
(which  we  will  caU  "e”  or  “ecoentricity”).  Once  these  perturbations  are  added,  the 
periodic  nature  of  the  orbit  may  no  longer  exist  even  in  theory.  It  is  hjrpothesised, 
however,  that  motion  very  near  the  three-body  orbit  will  still  be  at  least  quasi- 
pexiodic  and  that  the  appBcation  of  modal  control  will  stalnlise  the  satellite. 

S.l  The  Perturbed  System 

Let  us  return  to  the  Hamiltonian  for  the  system.  Without  rewriting  it,  we  note 
that  it  can  be  expanded  as  a  Taylor’s  series  about  the  periodic  orbit  (»e,  X  =  XpOt 
Sx  =  0).  Note  that  this  is  not  a  McLaurin  series  even  though  8z  —  0. 

H{X)  =  Ho  +  +  0(«x*)  (5.1) 


where 


Ho  =  H{Xp) 


The  first  term.  Ho,  is  a  scalar  constant  because  it  is  the  value  of  the  Hamiltonian 
for  the  periodic  orbit;  H\  is  identically  sero  because  it  describes  the  motion  of  the 
periodic  trajectory  with  respect  to  itself;  Hi  i>  the  quadratic  term  and  is  the  portion 
which  provided  the  unperturbed  problem  discussed  in  this  thesis  so  far  (the  Floqnet 
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pxoUem,  which  is  a  Uiiear  system  in  the  modal  vaiiaUM).  ^em  are  pertaihiag 
forces  out  there,  and  in  ord»  to  find  their  effects  on  the  system,  we  mnst  bcgia 
aipun  with  the  statement  of  the  new  problem.  When  we  have  stated  it  (see  section 
5.2),  and  extracted  its  mathematics,  it  wiU  be  seen  that  the  perturbations  can  be 
represented  by  a  Hamiltonian  of  their  own  (J7^t)  and  that  we  can  add  their  effects 
to  the  system  by  adding  in  the  perturbing  Hamiltonian. 

H  =  Ho  +  Hx  +  Hi  +  0{Sx*)  +  (5.3) 

Recall  from  Eqs.  2.29  and  2.34  that  we  have  already  accounted  for  the  second  order 
term  in  Eq.  5.1. 

Si  =  A{t)Sx  =  Z^Sx  (5.4) 

Oox 

We  made  the  assumption  that  that  the  third  order  term  was  smaller  than  the  per¬ 
turbations  so  we  solved  the  following  problem. 

«  =  +  (5.5) 

It  may  not  be  clear  yet,  but  this  is  a  first  order  linear  ODE  with  a  forcing  function. 
We  would  like  to  write  the  problem  in  the  modal  variables  for  simplification  so  we 
take  the  transformation  to  modal  variables  from  section  3.1. 

Sx{t)  =  (5.6) 


Taking  the  time  derivative  yields 

(5-7) 
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Now  combiaiiig  Eq/i.  5.4,  5.5  aad  5.7  we  5«ve 


+#•(<«(«)  =  ^‘)*»+z^  (M) 

Now  we  tabftiiete  {m  ike  modal  euiaklet  awd  aohre  for  a  diffenfttial  equlfoii  ia  f 
(let  as  also  drop  the  tiioe  iadex). 


ii^F-\AF-F}ii  +  F-^Z 


da 


(5.») 


RecalUag  Eq.  2.40 


F  =  AF-FJ 


(5.10) 


aad  solving  for  J 

J  =  F-\AF-F)  (5.11) 

we  see  we  can  simplify  Eq.  5.9  as  follows 

n  =  j,i+F-'z?^  (5.1*) 

It  is  now  easier  to  see  that  this  is  a  first  order  linear  ODE  with  a  ibrcing  innction. 

The  control  problem  that  was  solved  in  Chapter  IV  coaid  have  been  represented 
as  a  linear  system.  The  eqaations  of  motion  for  the  system  controlled  can  be  written 
in  the  modal  variables  as 

m = jm  (»ij) 

This  is  a  linear,  constant  coeflident  set  of  homogeneons  differential  eqaations  (re¬ 
member  that  this  system  resalted  from  aaafysis  of  a  linear  time-periodic  system  in 
the  state  variables).  In  the  pertarbed  system,  we  add  a  forcing  fhnction,  whidi  forces 
the  solation  to  be  the  snm  of  the  homogeneons  and  particnlar  solations.  The  system 
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(5.14) 


rnmlM  Kiir 


The  partkolar  loliitkm  hM  the  fotm  (5:  255) 

*(0  =  ."*(«,)  +  /  ,^‘->fl(r)4r  (5.16) 

but  we  do  not  need  the  aolntion.  We  need  a  ftatement  of  the  differential  equar 
tiona  to  allow  our  numerical  int^raUw  to  calculate  the  aolution  over  a  giTen  time 
internal.  In  order  to  write  the  differential  equation!,  we  need  an  expieaaion  fot 

/^jy _ 

B(t)  =  F'~^Z — The  next  aectkm  will  introduce  the  proUem. 
ox 

S.i  The  EUipUc  Reairieted  ProUem 

The  Reatricted  Problem  of  Three  Bodiea  haa  aeveral  modificationa.  Aa  Saebe- 
hely  put  it 

The  motion  of  the  primariea  mnat  aatiafy  the  differential  equation! 
goreming  the  dynamic!  of  two  gravitational  bodiea.  Conaeqnently  the 
jNnmariea  might  deacribe  dUptic,  parabolic  or  hyperbolic  orbita.  The 
apedal  caae  of  drcnlar  motion  ia  [the  restricted  three-body  problem].  Its 
aimplest  generalisation,  when  the  primariea  deacribe  elliptic  orbita  .  .  . 
ia  not  trivial  aince  its  conaequencea  are  [that  the  Hamihoniaa  will  be  a 
function  of  time](7:  587), 

It  ia  the  purpose  of  this  section  to  introduce  the  method  by  which  we  obtained  the 
Hamiltonian  and  Hamilton’s  equations  of  motion  ior  the  perturbations  introduced 
by  the  fact  that  the  Earth’s  orbit  has  aa  eccentricity  of  about  0.0167.  The  next 
section  will  discuss  the  added  effect  of  the  moon  which  orbits  the  Earth  with  a 
frequency  of  about  once  every  29.53  days  (sidereal).  Figure  5.1  shows  the  setup 
for  the  elliptic  problem  where  the  center  of  mass  symbol  on  the  left  indicates  the 
Earth>moon  barycenter.  The  distance  from  the  sun  to  the  Earth-moon  barycenter 
is  no  longer  constant,  and  is  represented  by  r  (not  to  be  confused  with  f).  Brouwer 
and  demence  (1:  76)  give  the  expansion  for  this  distance  in  small  eccentricity,  e.  To 
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I 


Figue  5.1.  The  Elfiptic  Beetricted  Problem 

first  order  it  is 

r  =  tt(l  —  ecos  jl#)  +  0(e*).  (6.16) 

We  take  a,  the  semi-major  aids,  to  be  1  AU,  the  same  as  befiote.  M  is  the  mean 
anomaly  {M  =  n(t  —  to))  where  the  mean  motion,  n,  is  eqnal  to  1,  so  that  Af  is 
eqniyalent  to  t.  Note  that  the  angular  velocity  of  the  rotating  reference  frame  is  no 

A 

longer  constant.  We  denote  w  =  vk.  Writing  to  first  order  from  the  expansion  for 
true  anomaly,  v  (1:  76), 

u  =  M  +  2e^M  +  0{e*)  (5.17) 

Tsking  the  time  derivative  yidds 

u;  =  i>  =  n-|-2encosAf -I- 0(e’)  (5.18) 
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Wb  Aotm  ooofdiaatet,  i,  f ,  i  tkal  aae  faded  by  tbe  bekveea  tbe  Eutb  aad 

the  ran. 

m  =  «/r 

y  =  f/r  (8*1®) 

J  =  xfr 

SolYuig  Ibr  tlie  spatial  coordiiiaies  and  mibatitiiting  ftcsti  Eq.  5.16  iot  tbe  'value  of 
yields 

c  =  £(1  — eoost) 

y  =  y(l  — eoost)  (8-20) 

X  =  i(l  — ecost) 

and  taking  derivatives  with  respect  to  time, 

X  =  2(1  —  e«)st)  +  enaint 

y  =:  y(l  —  ecost)  +  ey sin t  (8-21) 

i  =s  2(1  — ecost) +  eirint. 


Tbe  procedure  to  IbDow  when  writing  equations  in  scaled  coordinates  is  to  write 
in  spatial  coordinates  first  and  then  transCorm  into  tbe  scaled  coordinates 
itaing  equations  5.20  and  5.21.  Tbe  position  vector  is  as  before 


f  =  ari  +  yi  +  sb 


(8.22) 


in  rotating  spatial  (not  scaled)  coordinates.  Taking  tbe  inertial  time  derivative  gives 
tbe  inertial  velocity 

V  =  (2  —  y  —  2ey  cost)i  +  (y  +  »  +  2excost)j  +  ib  +  0{e^).  (8.23) 
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The  inertial  'vdocity  can  then  be  written  in  terms  ci  the  new  oocwdinates 


i(l  —  eco8t)  +  eisint  —  —  ecost)  —  2ey  cost) 

y  +  eysint  +  e(l  — ecost)  +  2eccost)  +C?(e*).  (5.24) 

i(l-ecosO  ^ 

Continuing  in  the  same  manner  as  developed  in  Glu4>ter  11  we  calculate  the  kinetic 
energy,  T. 

^  —  j((*  —  P  —  ey  cost  —  e*  cost  +  e*  sin  <)* 

+  (y +  *  +  c*cos<  — eycost  +  cpsint)*  (5.25) 

+  (i  -  ei  cos  I)*]  +  0(e*) 

&T 

The  generalised  momenta  can  be  obtained  from  pi  = 

oqi 

9T 

J>s  = -^r  ^  (*  — y  —  eycost  — c*cost  +  c*sint)(l  — ecost)  (5.26) 
dx 
OT 

Pg=:-^  «  (y -f  e  +  ex  cost  —  ey  cost -f- ey  sin  t)(l  -  e  cost)  (5.27) 

dy 


^‘^62 


(i  —  eicoat)(l  —  ecost) 


(6.28) 


which  are  not  the  components  of  the  inertial  velocity  vector.  Now  note  that  after 
some  algebra  the  kinetic  energy  can  also  be  written 


T=l 


2  [(1  — ecost)’  (1  — ecost)’  (1  — ecost)’J* 
This  relation  will  prove  useful.  It  will  also  be  useful  to  solve  for  the  qi's. 


(5.29) 


y  = 


_  _  pt  ,  y  +  cy  cos  t  —  e£  sin  t) 
*  (1  “  c  cos  t)’  (1  ~  c  cos  t) 

Pff  £  +  ex  cos  t  —  ey  sin  t) 


(1  —  ecost)’ 


(1  —  ecost)’ 


(1  —  ecost) 


(5.30) 
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Wlwa  it  is  zM&ed  that 


(1  —  ecost)”*  =  l-|-2ecost-i-C)(e’) 

(5.31) 

and 

(1  —  eoost)"'  =  1  -f>  ecost  +  0(e*) 

(5.32) 

then  the  qis  can  be  written 

i  =  Pi  +  y  +  2ey  cos  t  -f  2ep§  cos  t  —  ex  sin  t 

y  =  p^  —  X  +  2ex  cost -i- 2ep|  cost  —  ey  sin  t 

i  =  Pj  -1-  2epi  cos  t 

(5.33) 

The  Hamiltonian  for  this  problem  can  now  be  written  even  though  the  poten¬ 
tials  have  not  yet  been  calculated. 

H  = 

+ 

+ 

Krf + + ri) + »»  -  Jf*  +  Kp 
epjcosi  +  ej^cost  epjcosf 

2ep«ycost  —  2^xcost  —  ep§xcost  —  ep^ycost  + 

(5.34) 

The  first  line  of  Eq.  5.34  is  the  same  as  in  Chapter  II.  It  is  the  "restricted  problem” 
Hamiltonian,  Hrp.  The  rest  of  the  Hamiltonian  is  dne  entirely  to  the  perturbations. 
We  will  write  Hamilton’s  equations  of  motion  after  the  potentials  have  been  found. 

5.3  The  Elliptic  Four  Body  Problem 

Imagine  now  that  the  moon  orbits  the  Earth  in  a  circular  orbit  in  the  plane 
of  the  ecliptic  and  that  their  barycenter  orbits  the  sun  in  a  slightly  eccentric  orbit. 
A  graphic  representation  of  the  problem  can  be  seen  in  Figure  5.2.  This  introduces 
another  perturbation  arising  from  the  moon’s  monthly  drcnit  of  the  Earth.  The 
potentiab  are  affected  by  this  motion  of  the  moon,  but  the  kinetic  energy  part  of 
the  Hamiltonian  is  not.  In  order  to  write  the  potentials  for  this  problem,  we  draw 
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Figure  5.2.  Frame  of  Reference  for  the  Restricted  Four  Body  Problem 

another  picture.  Figure  5.3  shows  the  vectors  to  the  satellite  from  each  of  the  three 
primary  bodies,  as  well  as  some  other  pertinent  distances.  The  angular  frequency  of 
the  moon  is  W0.  The  potential  for  the  satellite  is  written 

V  =  Vq  +  V.  +  V^ 

_  (5.35) 

»‘0.s  *0,8 

where  the  subscripts  denote  the  body  from  which  the  potential  originates:  O  is  the 
sun,  0  is  the  Earth,  and  0  is  the  moon.  The  masses  need  to  be  parameterised 
to  keep  the  units  consistent,  so  we  define  several  |»'s  to  represent  them.  These  are 
not  to  be  confused  with  the  original  /i  for  the  three-body  problem,  though  they  are 

3.00348058  x  lO"* 

3.69429133  x  lO"*  (5.36) 

1.21505814  X  10-> 


similar. 


Me 

Ma 


m® 


n»©  +  w»®  + 
na®  +  w»®  +  wa® 

TO®  -hna® 
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Figure  5.3.  More  Detail  on  the  Beatricted  Four  Body  Problem 

The  quantity  /ij  is  the  parameter  for  the  three  body  problem  involving  the  Earth 
and  moon  as  primaries. 

In  order  to  separate  the  perturbed  Hamiltonian  from  the  rest  of  the  Hamilto¬ 
nian,  it  will  be  necessary  to  expand  the  distances  in  the  denominators  of  the  potential 
in  the  small  parameters  (eccentricity  of  Earth’s  orbit,  e,  and  the  semimajor  axis  of 
the  moon’s  orbit,  a^).  We  will  follow  this  process  closely  for  the  first  term  (the  sun’s 
potential)  and  the  other  terms  (Earth’s  and  moon’s  potentials)  will  be  derived  in 
less  detail.  We  are  careful  to  write  the  vectors  in  spatial  coordinates  first  and  then 
transform  them  into  scaled  coordinates. 

^©,3  =  (*-#*»■)»  + »i  +  **  (5.37) 

=  (*  -  (1  -  e  cost)^]  i  +  yj  -I-  xh  +  0(e*)  (5.38) 
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We  are  expanding  iti  magnitude  in  the  denominator,  lo  we  write 

=  'ii  =[(•-/•  +  (««).«)»+»>  +  .>]■*  +  0(e*)  (6J«) 

Now  we  convert  to  scaled  coordinates. 


(S.«) 

Simplifying,  and  gronping  the  large  terms  from  the  small  ones  gives 


’•ils  =  [(*-/*)*  +  »*  +  »’ 

—  2e(x  —  /i)*  cos  t  —  2ep’  cost  —  2ei’  cost]~i  +  0(e*) 
Now  we  expand  to  first  order  using  a  binomial  expansion  to  obtain 


••si  =  [(*-/.)• +»•+?)•* 

X  2e  cos  <(-(*  —  /»)*  --  y*  -  **)  +  0{e^) 


(6.41) 


(6.42) 


which  can  be  simplified  greatly  by  letting  pi  =  ^{x  —  /i)*  +  y*  +  i*  (as  in  Eq.  2.6) 
so  that 


’’©i  —  + 


Pi 


ecost 


Pi 


(6.43) 


Now  we  have  an  expression  for  the  first  term  of  Eq.  6.36;  noting  that  the  mass  of 
the  sun  is  (1  —  p)  we  have 


(1-m)  {l-p)ecoit 
Pi  Pi 


(6.44) 

(6.46) 


where  the  subscript  rp  denotes  the  restricted  problem  (RTBP)  and  pert  denotes 
"perturbed”  so  that  the  potential  is  the  sum  of  the  potential  due  to  the  three  body 
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problem  and  tbe  potential  dne  to  pertnrbatkas.  F<»  the  ran  potential,  ecoentncity 
is  the  only  perturbation,  so  does  not  appear. 

The  potentials  due  to  the  Earth  and  ino<m  can  be  scdved  in  a  similar  way.  The 
vectors  to  the  satellite  from  the  Earth  and  moon  in  spatial  coordinates  are 

^•,8  =  [*  +  (1  -  /*)(!  -  e  cos  I)  +  a^fii  cos((«;0t)]  i 

+  [y  +  «c/*a  «in(«0Oli  +  ** 

(5.46) 

^04  =  (*  +  (l-M)(l-eco«O-«0(l-Ma)«>*("0O]» 

+  [y  -  «0(1  -  /*a)«“("0O]i  +  ** 

So  the  inverse  distances  are 

’“ei  =  {{*  +  (l-/i)(l-cco8t)  +  O0/tacos(w0<)]* 

+  [y  +  «0/*3  «in(w0t)]’  +  z* 

(5.47) 

’•ojs  =  {l*  +  (l-M)(l-eco8<)-®0(l-M»)«»(‘^0O]* 

+  [y  ~  “0(1  -  Ma)  8in(«0t)]’  +  z*  }“ » 

We  then  transform  into  scaled  coordinates  from  Eq.  5.20  to  give 

»’©,s  =  {[*(1 -c cost) +  (1 -/*)(! -c cost) +  O0;»acos(<i;0t)]* 

+  [jf(l  -  ccost)  +  00/1]  sin(fa;0t)]’  +  [z(l  -  ccost)]*  }"^ 

(5.48) 

’•0?3  =  {[i(l  -  e cos t) +  (1- /*)(!- ccost) -00(1- Pa) cos(u;0t)]* 

+  [y(l  -  ccost)  -  00(1  -  pa)8bi(w0Of  +  “  ecost)]*  }“» 

These  are  expanded  in  both  c  and  00,  and  the  second  order  terms  are  removed. 
We  again  let  pa  =  y/{x  +  1  —  p)’  +  y*  +  z*  and  multiply  by  the  respective  masses  to 
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obtain 


Vi  = 


Vo  = 


-£*  +  ^  =  v*,+v, 

pi  n 

I  _  y  ,  y 

^  ^  -Ve,  +  V, 


(5.49) 


where  the  /3*s  axe  given  by 


=  --/>acoete  +  /t3[y«in(wet)  +  (£  +  l-/t)coa(<«'e01^ 
fia  =  coate +  + 


(5.50) 


Again  we  see  that  the  potentials  are  the  stun  of  the  previous  potential  and  a  new 
perturbation  potential. 

In  order  to  get  equations  of  motion  for  the  system,  we  take  the  derivative  of 
the  Hamiltonian  with  respect  to  the  state  vector  as  described  in  Chapter  11 


X  =  Z 


dX 


(5.51) 


Rather  than  repeat  the  work  listed  in  Chapter  II,  we  will  give  the  partial  derivatives 
of  only  the  perturbed  Hamiltonian,  Hp,  with  respect  to  the  elements  of  the  state 
vector,  where  we  will  let  the  subscript  p  dmiote  perturbed. 


Hp  =  eplcoBt  + ep- cost +  ep^ccmt  +  2epty  cost  — 2epgx  cost 
-  epsxidat- epgy^t  +  Vop  +  V^  +  V^p 


(5.52) 
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dH, 

dVx, 

^dv^. 

dx 

II 

1 

1 

s 

1 

e 

«« 

dx 

^  dx 

dH, 

dy 

=  2epi  cos  t  —  eyy  sin  i  + 

.  dV^ 

dy  ^  dy  ^ 

dv^ 

dy 

dH, 

dz 

dv^  av«p  av0, 

dz  dz  dz 

dH, 

dp» 

=  2ey  cost +  2ep«  cost  — 

exsint 

dH, 

dpfj 

=  — 2ex  cos  t  +  2e|iy  cos  t 
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The  partials  of  the  potential  terms  ate  given  by 
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dH,  . 


We  now  have  onr  expression  for  in  the  variational  state  variables.  If  we 

dH 

premnltiply  it  by  Z  and  add  it  to  Z-^  we  get  the  right  hand  side  of  the  total 

in 

problem  (see  Eq.  5.5).  A  premnltiplication  by  F~^  will  transform  the  problem  into 
the  modal  variables  (see  Eq.  5.12). 


The  total  problem  has  been  stated  and  the  integrator  will  integrate  these  equa¬ 
tions  of  motion  readily.  It  should  be  noted  here  that  the  introduction  of  a  forcing 
function  to  onr  linear  system  creates  a  particular  s<^ution,  which  is  added  to  the  ho¬ 
mogeneous  solution  to  form  the  total.  The  total  solution  is  the  one  that  is  observed, 
and  it  is  the  one  whose  modes  are  controlled.  The  nature  of  the  particular  solution, 
which  may  have  some  significance  in  the  control  of  the  satdlite,  is  ignored.  It  is 
known  that  if  the  unstable  mode  of  the  particular  solution  diverges  in  one  direction 
(say  toward  +oo)  then  the  homogeneous  solution  can  be  made  to  diverge  in  the 
opposite  direction  (a  very  large  negative  number)  so  that  their  sum  is  quite  small. 
The  consequences  of  this  simultaneous  divergence  during  the  control  of  the  satellite 
have  not  yet  been  discovered. 
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VI.  Reaulta  of  Perturbation  Controller 

The  problems  associated  with  the  perturbed  orbit  were  many.  Tbe  periodic 
orbit  no  longer  existed  in  the  same  sense  that  it  did  behne.  The  satdfite  placed  near 
the  three-body  orbit  was  able  to  nearly  ckwe  on  itself  several  times  be&ne  diverging 
(see  section  4.1),  but  the  perturbations  cause  the  satdUte  to  diverge  much  socmer  — 
before  completing  a  single  orlnt.  See  Figure  6.1. 

6.1  SimpUfication  of  System 

When  perturbations  from  the  eccmitricity  and  the  moon  were  included  in  the 
model,  the  controller  failed  to  adequately  control  the  sateUite  in  the  orbit.  The 
reasons  for  this  are  not  precisely  known.  Perhaps  the  second  order  effects  should  be 
taken  into  account  before  the  mode  shapes  will  be  similar  enough  to  the  unperturbed 
case  for  the  controller  to  be  effective.  Perhaps  the  difference  betsreen  the  particular 
and  the  homogeneous  solutions  caused  the  total  solution  to  compound  numerical 
errors  so  that  the  simulation  broke  down.  In  any  case,  the  pe^nrbations  due  to 
eccentricity  were  put  aside,  and  the  effects  of  the  moon  alone  were  studied.  The 
results  shown  in  this  section  will  be  the  results  of  the  controller  acting  upon  a 
simulated  system  where  the  earth  is  in  a  drculai  orbit  about  the  sun,  but  the  moon 
is  also  in  orbit  about  the  earth,  and  is  affecting  the  satellite  with  its  gravity.  The 
equations  of  motion  and  the  computer  code  which  integrates  them  remain  exactly 
the  same,  with  the  simple  change  of  setting  e  equal  to  zero. 

With  the  simplified  system,  the  modal  control  scheme  was  successful.  Even 
with  e  =  0,  the  uncontrolled  satellite  diverged  before  completing  a  single  orbit,  but 
the  divergence  was  somewhat  slower.  See  Figure  6.2.  Using  modal  control,  it  was 
kqpt  in  the  vicinity  of  the  orbit  indefinitely  as  long  as  the  threshold  was  small  enough. 
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Figoie  6.3.  Average  Station-Keeping  Coat  for  Moon-Pertnrbed  Orbit,  Threshcdd  = 
10-* 


6.S  Rtsuha 

Let  U8  look  firat  at  a  threshold  which  gives  a  very  similar  result  to  the  un¬ 
perturbed  case.  The  threshold  in  this  case  was  about  10~*.  See  figures  6.3,  6.4 
and  6.5.  All  the  simulations  in  this  section  were  run  for  31  time  units  (canonical) 
which  is  close  to  10  orbital  periods  or  about  5  years.  The  station-keeping  cost  for 
this  first  simulation  is  very  high.  When  the  canonical  units  ate  converted  to  more 
conventional  units,  it  cost  about  200  m/sec  per  year  to  control  the  satellite  with  this 
threshold. 
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Figure  6.4.  Plot  of  the  Real  Modes,  Threshold  =  10  * 

The  mode  shapes  for  this  simulation  were  somewhat  regular  (figure  6.4).  Even 
though  the  maneuvers  do  not  sero  out  both  modes  (the  stable  mode  is  given  a 
positive  value  while  the  unstable  mode  is  set  to  near  sero-just  like  in  section  4.4), 
the  behavior  seems  to  be  fairly  predictable.  There  are  some  funny  curves  in  the  plot, 
and  these  are  due  to  the  perturbations. 

Figure  6.5  shows  the  approximate  ‘‘Earth  view”  for  the  orbit.  It  is  a  trace  of 
the  Y  and  Z  positions  of  the  satellite  throughout  the  simulation.  It  shows  that  there 
is  considerable  wobble  in  the  orbit,  which  may  or  may  not  be  acceptable  depending 


Figure  6.5.  Plot  of  the  Y  -  Z  Projection  of  the  Orbit,  Threshold  =  10“* 

on  the  size  of  the  sun’s  radio  disk  and  other  constraints.  Interestingly  enough,  the 
orbit  is  both  cheaper  to  control  and  more  controlled  at  this  threshold  than  the 
unperturbed  orbit  was.  Compare  to  Figure  4.10. 

If  we  decrease  the  threshold  for  the  unstable  mode,  we  can  decrease  the  control 
cost  somewhat.  Figure  6.6  shows  the  cost  for  a  threshold  of  10~*.  Although  the 
total  average  cost  is  lower  for  this  option  (about  20  m/sec  per  year),  it  requires  very 
frequent  maneuvers.  During  this  simulation,  over  300  maneuvers  were  performed, 
averaging  about  5  per  month. 
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Time  (Csmonicfd  Units) 


Figure  6.6.  Station-Keeping  Cost  with  Threshold  =  10~*  Showing  the  Maneuvering 
is  Extremely  Frequent. 
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Figure  6.7.  Y  -  Z  Projection  (Like  the  View  from  Earth)  of  Controlled  SatdUte’s 
Orbit,  Threshold  =  10~* 

The  Y  -  Z  projection  of  the  orbit  (see  Figure  6.7)  shows  considerable  tightness 
in  the  control  of  the  orbit.  The  slight  anomalies  that  can  be  seen  hme  are  due  to 
the  periodic  fluctuations  in  the  perturbing  force  of  the  moon’s  gravity  and  not  to 
divergence  and  subsequent  control  of  the  unstable  mode. 

Since  the  manuevers  were  so  frequent  with  the  threshold  at  10~*,  an  attempt 
was  made  using  trial  and  error  (shall  we  say,  “trial  and  success”)  to  optimise  the 
threshold  to  decrease  the  number  of  maneuvers  required  and  to  decrease  the  total 
station-keeping  cost  if  possible.  The  attempt  was  successful  with  a  threshold  of 
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Figure  6.8.  Station-Keeping  Cost,  Threshold  =  2  x  10~' 

2  X  10~'‘.  Figure  6.8  shows  the  cost,  which  is  about  10  m/sec  per  year  when  units 
are  converted.  It  can  also  be  seen  that  there  are  fewer  maneuvers  —  about  one  per 
month.  The  average  maneuver  was  about  72.5  cm/sec. 

Figure  6.9  is  a  plot  of  the  real  modes.  It  seems  rather  jumbled,  but  this  is 
because  the  threshold  is  on  the  same  orda  (spatial  and  temporal  frequency)  as  the 
gyrations  of  the  modes  so  that  the  maneuver  occurs  in  unpredictable  positions  on 
the  plot.  What  is  important  is  that  the  unstable  mode  never  exceeds  the  threshold, 
and  that  it  is  allowed  to  “wander”  until  its  wandering  begins  to  diverge. 
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Stable  Mode 


Unstable  Mode 

Figure  6.9.  Real  Modes  Under  Control  for  Threshold  =  2  x 
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Figure  6.10.  Y  -  Z  Projection  for  the  Satellite  Path,  Threshold  =  2  x  10~* 


The  Earth-view  for  this  threshold  shows  a  very  small  divergence  from  the 
nominal  orbit,  which  would  be  acceptable  under  most  circumstances.  If  the  path 
crossed  too  dose  to  the  sun,  the  best  alternative  would  be  to  find  a  new  periodic 
reference  orbit  in  the  three-body  problem  which  had  larger  amplitude  rather  than 
attempt  tighter  control. 

Let  us  look  at  the  cost  as  a  function  of  the  threshold  for  the  perturbed  satellite. 
Figure  6.11  shows  the  station-keeping  cost  for  a  number  of  different  thresholds. 


Average  Station-Keeping  Cost  (cm/sec  per  Year) 


Threshold  Value 


Figure  6.11.  Station-Keeping  Cost  as  a  Function  of  Threshold 
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Figure  6.12.  Cost  Comparison  for  Perturbed  (+)  vs  Unperturbed  (o)  Systems 

It  is  easy  to  see  that  the  perturbed  satellite  is  more  expensive  to  control  than 
the  imperturbed  satellite  by  looking  at  Figure  6.12.  When  two  comparable  scenarios 
are  compared  (threshold  =  10~^  for  the  unperturbed  and  2  x  10'*^  for  the  perturbed), 
we  see  that  the  perturbed  satellite  requires  1000  times  more  control  energy  than  the 
unperturbed  satellite. 


6-13 


6.3  Behavior  of  the  Other  Modes 


Once  again,  the  behavior  of  the  otha  modea  nnder  control  ia  intereating  and 
will  therefore  be  included  here.  The  perturbing  effect  of  the  moon  cauaed  aome 
intereating  ahapea  in  them. 

Let  ua  look  at  Figure  6.13.  During  each  of  the  manenveta,  there  ia  a  diaconti- 
nuity  in  the  plot,  plotted  with  a  atraight  Hne,  but  hardly  viaible  on  thia  plot  becauae 
the  manuevera  have  auch  amall  effect  on  thia  mode.  Between  maneuvera,  when  the 
aatellite  ia  drifting  freely,  the  imaginary  modea  tend  to  trace  out  a  akewed  5-pointed 
atar  ahape,  which  alowly  drifta  away  from  the  ori|^.  Thia  ia  due  to  the  moon’a 
gravity,  and  aa  we  would  expect,  the  period  of  the  perturbation  effecta  ia  about  aix 
timer  the  period  of  the  orbit.  Thia  maker  aenae  aince  the  moon’a  motion  haa  a  period 
of  one  month  to  the  halo  orbit’a  aix  montha. 

Figure  6.14  ia  like  the  previoua  figure,  with  the  aero  modea  plotted  inatead  of 
the  imaginary  modea.  When  the  aateUite  ia  drifting  freely,  the  value  of  the  modea 
ahowa  the  atar  ahape,  and  when  a  maneuver  ia  performed,  their  value  changer  alightly 
in  aome  predictable  direction,  to  repeat  the  figure  until  the  next  maneuver. 

6.4  Summary 

The  aatellite  perturbed  by  the  eccentricity  war  not  controllable  with  thia  con¬ 
troller.  When  the  moon  only  war  included,  however,  the  aituation  improved  dra¬ 
matically.  The  modea  diaplayed  their  expected  behavior  within  reaaon,  with  aome 
unpredictable  perturbation  effecta.  The  aatellite  war  controllable  with  reaaonable,  if 
higher,  control  coata.  The  reaaona  for  the  failure  with  the  eccentricity  included  are 
diacuaaed  in  Chapter  VII. 
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Mode  Four 


Mode  Six 


Figure  6.14.  Zero  Modes  Plotted  Together  While  Sfttellite  is  Under  Realistic 
Control 


VII.  Conelusioru  and  RecommendaHons 

7.t  Conduaiona 

Modal  control  works  wdl  on  the  unperturbed  three-body  problem  satdhte. 
Using  maneuvers  calculated  by  observation  of  the  modes,  the  satellite  was 

contrdled  for  as  long  as  desired.  The  concept  was  simple,  and  the  ^pHcatkin  was 
ideal.  Modal  control  held  promise  of  being  more  efficient  than  other  c^timal  control 
algorithms  for  this  type  of  orbit. 

We  learned  from  our  cheater  controller  that  each  of  the  modes  had  a  charac¬ 
ter  which  depended  on  its  Poincare  exponent.  The  character  was  exactly  what  was 
expected  as  long  as  the  magnitude  remained  small  (within  the  Unaar  region).  The 
real  exponents  corresponded  to  an  exponentially  decreasing  mode  (negative  expo¬ 
nent)  and  an  exponentially  increasing  mode  (positive  teal  exponent).  The  imaginary 
modes  (corresponding  to  the  imaginary  exponents)  oscillated  sinusoidally  and  were 
90  degrees  out  of  phase  with  each  other  as  indicated  by  Figure  4.2,  where  they  form 
a  circle  when  plotted  together.  The  aero  modes  were  static  to  within  reasonable 
ntimerical  accuracy. 

The  efficiency  of  the  controller  on  the  unperturbed  satellite  was  improved  by 
decreasing  the  threshold  value  up  to  a  point.  After  that,  the  cost  actually  increased 
when  the  threshold  was  decreased.  This  was  due  to  the  numerical  inaccuracy  of  the 
integration  at  such  small  values  of  the  modes  (less  than  10~*).  The  station-keeping 
costs  for  the  unperturbed  satellite  were  phenomenally  low.  For  a  15  year  mission, 
the  satellite  could  be  maintained  with  1.5  cm/sec  of  AV  total  (threshold  =  10~^) 
with  2  small  maneuvers  per  year. 

The  perturbed  satellite  had  a  different  character.  The  perturbation  due  to 
eccentricity  was  very  large,  and  very  probably  pulled  the  satellite  out  of  the  linear 
region  very  quickly,  which  invalidated  the  shapes  of  the  modes  and  thus  the  con¬ 
trol  scheme.  The  perturbation  due  to  the  moon  was  not  as  large,  and  the  modes. 
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wliile  noticeably  different  than  the  unperturbed  modes,  were  still  recogniaeaUe.  We 
were  not  aUe  to  decrease  the  station-keeping  cost  nearly  as  much  by  decreasing  the 
threshold.  This  makes  sense,  since  decreasing  the  threshcdd  caused  a  lot  of  unnec- 
cessary  maneuyering.  Tim  moon  was  very  quick  to  pull  the  satellite  off  the  periodic 
orbit  and  thus  to  cause  the  modes  to  grow,  but  only  to  a  point,  and  if  we  maneuyered 
right  away,  we  missed  the  fact  that  the  moon  just  as  sutdy  “pushed”  the  satdhte 
tight  back  near  the  orbit  during  the  second  half  of  the  month.  The  best  strategy 
was  to  wait  to  maneuyer  when  the  unstable  ^aode  was  actually  beginning  to  grow 
of  its  own  accord,  and  not  due  the  perturbation.  The  threshold  that  worked  best 
for  the  perturbed  satellite  was  2  x  10~*.  The  station-keeping  cost  for  this  threshold 
was  approximately  1  m/sec  per  year,  nearly  3  orders  of  magnitude  higher  than  the 
unperturbed  cost. 

7.2  Recommendaiiona 

The  perturbations  due  to  eccentricity  must  be  included  in  any  model  that 
represents  itself  as  realistic.  There  are  other  perturbations,  but  they  are  snudler 
(2).  It  is  recommended  that  some  scheme  for  including  this  perturbation  be  found, 
which  allows  the  mode  shapes  to  remain  apparent,  so  that  modal  control  will  work. 
Perhaps  a  model  based  on  an  expansion  to  second  order  in  the  eccentricity  would 
work. 

In  order  to  yalidate  the  modal  control  scheme  as  workable,  there  must  also  be 
studies  done  on  the  probabilistic  problem  of  orbit  (and  thus  mode)  determination. 
The  halo  satellite  presents  special  problems  in  this  area,  which  haye  been  studied  by 
Howell  and  Pemicka  —  see  references  (3)  and  (4),  but  not  with  modal  control. 

The  consequences  of  the  particular  solution  diyergence  while  the  total  solution 
remains  finite  haye  not  been  addressed.  If  this  scheme  is  to  be  workable,  the  solutions 
to  the  mathematical  model  must  be  yery  well  understood,  so  it  is  recommended  that 
some  effort  be  put  forward  in  this  endeayor. 
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were  low,  and  the  significance  is  that  the  controller  did  nothing  to  suppress 
modes  that  were  oscillatory  or  converging.  The  effect  of  the  moon's  motion  allowed 
the  scheme  to  operate  with  reasonable  control  costs,  but  the  effect  of  eccentricity 
caused  divergence  in  spite  of  the  controller. 
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